{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping and Plotting the ICESat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "import h5py\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/icesat/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.H5'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### process it & subsetting\n",
    "def glah14_to_df(filename,bounds):\n",
    "    ## sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "    ## Bounds are [Longitude_West, Longitude_East, Latitude_South, Latitude_North]\n",
    "    f = h5py.File(filename, 'r')\n",
    "    lat = f['Data_40HZ']['Geolocation']['d_lat'][:]\n",
    "    lon = f['Data_40HZ']['Geolocation']['d_lon'][:]\n",
    "    elev = f['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]\n",
    "    sec = f['Data_40HZ']['Elevation_Corrections']['d_satElevCorr'][:]\n",
    "    scf = f['Data_40HZ']['Quality']['sat_corr_flg'][:]\n",
    "    satndx = f['Data_40HZ']['Quality']['i_satNdx'][:]\n",
    "    dem = f['Data_40HZ']['Geophysical']['d_DEM_elv'][:]\n",
    "    \n",
    "    glah14_df = pd.DataFrame({'Latitude':lat,'Longitude':lon,'Elevation':elev,\n",
    "                            's_El_Corr':sec, 's_Corr_f':scf,'in_sat':satndx,\n",
    "                            'DEM':dem})\n",
    "    #### Subsetting\n",
    "    glah14_df_subset = glah14_df.loc[(glah14_df['Longitude']>=bounds[0]) \n",
    "                          & (glah14_df['Longitude']<=bounds[1])\n",
    "                          & (glah14_df['Latitude']>=bounds[2])\n",
    "                          & (glah14_df['Latitude']<=bounds[3])]\n",
    "    return glah14_df_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Spatial Bounds: \n",
    "tsl_sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "### test track of ICESat\n",
    "test_df = glah14_to_df(files[1],tsl_sp_ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         Latitude   Longitude  Elevation  s_El_Corr  s_Corr_f  in_sat    DEM\n",
      "159517  13.286178  104.542759     59.483      0.000         0       0  60.30\n",
      "159518  13.284618  104.542540     56.767      0.000         0       0  60.43\n",
      "159519  13.283058  104.542322     58.253      0.000         0       0  60.55\n",
      "159520  13.281498  104.542106     54.924      0.000         0       0  60.65\n",
      "159521  13.279939  104.541892     48.373      0.000         0       0  60.75\n",
      "...           ...         ...        ...        ...       ...     ...    ...\n",
      "160097  12.382718  104.416207      7.871      0.173         2       3  11.87\n",
      "160098  12.381162  104.415988      8.851      0.115         2       2  12.18\n",
      "160099  12.379608  104.415769      7.519      0.000         0       0  12.51\n",
      "160100  12.378055  104.415549      6.656      0.996         2      20  12.67\n",
      "160101  12.376502  104.415327      9.499      0.000         0       0  12.76\n",
      "\n",
      "[585 rows x 7 columns]\n"
     ]
    }
   ],
   "source": [
    "print(test_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ff2251962fcf49a491d3ba36c8b1a57b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Tracks on TSL')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(6,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.set_extent(tsl_sp_ex, crs=ccrs.PlateCarree())\n",
    "plt.scatter(test_df['Longitude'], test_df['Latitude'],s=1)\n",
    "\n",
    "request = cimgt.Stamen('terrain-background')\n",
    "ax.add_image(request, 9)\n",
    "plt.title(\"Tracks on TSL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "504b9c0315584103983adf258915dbe4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(test_df['Latitude'],test_df['Elevation'],'.',markersize=0.25, label='all heights')\n",
    "h_leg=ax.legend()\n",
    "plt.title('Heights from ICESat')\n",
    "ax.set_xlabel('Latitude')\n",
    "ax.set_ylabel('Elevation, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mapping the track with elevation information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bcc36cc087914b098994fddde83ef154",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(5,5))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.set_extent(tsl_sp_ex, crs=ccrs.PlateCarree())\n",
    "plt.scatter(test_df['Longitude'], test_df['Latitude'], s=2, c=test_df['Elevation'], alpha=.7, transform=ccrs.PlateCarree(), cmap='terrain')\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation (m)')\n",
    "\n",
    "##load Google Sat.Map\n",
    "#request = cimgt.GoogleTiles(style='satellite')\n",
    "request = cimgt.Stamen('terrain-background')\n",
    "ax.add_image(request, 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Removing outliers using flag information in itself"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
